First I need to load up the packages I’ll need
library(sf)
library(ggplot2) #development version!
## devtools::install_github("tidyverse/ggplot2")
library(tidyverse)
library(readr)
library(cowplot)
library(sp)
library(gridExtra)
library(dplyr)
library(ggrepel)
library(plyr)
Now I import my data. I filter for the Arran postcodes, (since Arran all begins ‘KA27’).
#Add download commands for data.
## Finding the Arran coordinates
arrancoordinates <- read.csv("../alldata/ukpostcodes.csv") %>%
filter(substr(postcode,1,4)=="KA27")
#Find way to replace with existing SIMD shape files
arransubsect <- read_sf("../alldata/Scotland_pcs_2011") %>%
filter(substr(label,1,4)=="KA27")
Now I load the SIMD data, containing the geometries (shapefiles) and SIMD data (percentiles, etc)
reorderedvector<- c("S01011174", "S01011171", "S01011177", "S01011176", "S01011175", "S01011173", "S01011172" )
arran2016 <- read_sf("../alldata/SG_SIMD_2016")[c(4672,4666,4669,4671,4667,4668,4670),] %>%
slice(match(reorderedvector, DataZone))
Arrandz2012 <- c(4409,4372,4353,4352,4351,4350,4349)
arran2012 <- read_sf("../alldata/SG_SIMD_2012")[Arrandz2012,]
arran2009 <- read_sf("../alldata/SG_SIMD_2009")[Arrandz2012,]
arran2006 <- read_sf("../alldata/SG_SIMD_2006")[Arrandz2012,]
arran2004 <- read_sf("../alldata/SG_SIMD_2004")[Arrandz2012,]
sharedvariables <- intersect(colnames(arran2016), colnames(arran2012)) %>%
intersect(colnames(arran2009)) %>%
intersect(colnames(arran2006)) %>%
intersect(colnames(arran2004))
arran20162 <- arran2016 %>%
select(sharedvariables) %>%
mutate(year="2016")
arran20122 <- arran2012 %>%
select(sharedvariables) %>%
mutate(year="2012")
arran20092 <- arran2009 %>%
select(sharedvariables) %>%
mutate(year="2009")
arran20062 <- arran2006 %>%
select(sharedvariables) %>%
mutate(year="2006")
arran20042 <- arran2004 %>%
select(sharedvariables) %>%
mutate(year="2004")
arransimd <- rbind(arran20162,arran20122,arran20092,arran20062,arran20042) %>%
mutate(
lon = map_dbl(geometry, ~st_centroid(.x)[[1]]),
lat = map_dbl(geometry, ~st_centroid(.x)[[2]])
)
arransimd$listID <- revalue(arransimd$DataZone,
c("S01004409"="S01004409/S01011174", "S01004372"="S01004372/S01011171", "S01004353"="S01004353/S01011177", "S01004352"="S01004352/S01011176", "S01004351"="S01004351/S01011175", "S01004350"="S01004350/S01011173", "S01004349"="S01004349/S01011172", "S01011174"="S01004409/S01011174", "S01011171"="S01004372/S01011171", "S01011177"="S01004353/S01011177", "S01011176"="S01004352/S01011176", "S01011175"="S01004351/S01011175", "S01011173"="S01004350/S01011173", "S01011172"="S01004349/S01011172"))
Now I want to overlay the postcodes by Datazone. To do this I’ve converted both the Arran coordinates and Arran (2016) shapefiles into spatial points/polygons, converted them into a common CRS, and then compared them by using ‘plyr::over()’. This gives me the object ‘namingdzpostcode’, with the postcode rows grouped into IDs (unidentified datazones).
simple.sf <- st_as_sf(arrancoordinates, coords=c('longitude','latitude'))
st_crs(simple.sf) <- 4326
exampleshapes <- sf:::as_Spatial(arran2016$geometry) %>%
spTransform(CRS("+proj=longlat +datum=WGS84"))
examplepoints <- sf:::as_Spatial(simple.sf$geom) %>%
spTransform(CRS("+proj=longlat +datum=WGS84"))
namingdzpostcode <- over(exampleshapes, examplepoints, returnList = TRUE)
I can then take a member reference from the orginal postcode list, which gives me a selection of the rows in that DZ. For simplicity I’ve written this as a new function. ##Mutate arrancoordinates to label the IDs
function100 <- function(argument)
{
argument <- arrancoordinates[namingdzpostcode[[argument]],] %>% mutate(DataZone=argument)
}
arrancoordinates <- lapply(1:7,function100)
arrancoordinates <- rbind(arrancoordinates[[1]], arrancoordinates[[2]], arrancoordinates[[3]], arrancoordinates[[4]], arrancoordinates[[5]], arrancoordinates[[6]], arrancoordinates[[7]])
arrancoordinates$listID <- revalue(as.character(arrancoordinates$DataZone),
c('1'="S01004409/S01011174", '2'="S01004372/S01011171", '3'="S01004353/S01011177", '4'="S01004352/S01011176", '5'="S01004351/S01011175", '6'="S01004350/S01011173", '7'="S01004349/S01011172"))
names(namingdzpostcode) <- c(unique(arransimd$listID))
///
library(rgdal)
library(leaflet)
library(ggmap)
postcodelist <- paste(unique(arrancoordinates$listID), "Postcodes", sep=" ")
datazonelist <- paste(unique(arrancoordinates$listID), "Datazones", sep=" ")
m = leaflet() %>% addTiles() %>% setView(-5.227680, 55.582338, zoom = 10)
Inputing example markers.
cliniccoordinates <- read.csv("../alldata/clinics.csv") %>%
dplyr::left_join(arrancoordinates, by="postcode")
Column `postcode` joining factors with different levels, coercing to character vector
#change to character
cliniccoordinates$X <- as.character(cliniccoordinates$X)
exampleshapes2 <- as(arransimd, "Spatial") %>%
spTransform(CRS("+proj=longlat +datum=WGS84"))
functionmap <- function(argument, argument2)
{
pal2 <- colorNumeric(palette = "viridis",
domain = argument)
listlistlist <- paste(datazonelist, argument, sep=" ")
m %>%
#alldatazones
addPolygons(data=exampleshapes2[exampleshapes2$year == 2004, ],
weight = 2,
label = listlistlist[29:35],
group = "2004",
fillOpacity =0.8,
color = ~pal2(argument[29:35]),
highlightOptions = highlightOptions(color = "black", weight = 2,
bringToFront = TRUE)) %>%
hideGroup("2004") %>%
addPolygons(data=exampleshapes2[exampleshapes2$year == 2006, ],
weight = 2,
label = listlistlist[22:28],
group = "2006",
fillOpacity =0.8,
color = ~pal2(argument[22:28]),
highlightOptions = highlightOptions(color = "black", weight = 2,
bringToFront = TRUE)) %>%
hideGroup("2006") %>%
addPolygons(data=exampleshapes2[exampleshapes2$year == 2009, ],
weight = 2,
label = listlistlist[15:21],
group = "2009",
fillOpacity =0.8,
color = ~pal2(argument[15:21]),
highlightOptions = highlightOptions(color = "black", weight = 2,
bringToFront = TRUE)) %>%
hideGroup("2009") %>%
addPolygons(data=exampleshapes2[exampleshapes2$year == 2012, ],
weight = 2,
label = listlistlist[8:14],
group = "2012",
fillOpacity =0.8,
color = ~pal2(argument[8:14]),
highlightOptions = highlightOptions(color = "black", weight = 2,
bringToFront = TRUE)) %>%
hideGroup("2012") %>%
addPolygons(data=exampleshapes2[exampleshapes2$year == 2016, ],
weight = 2,
label = listlistlist[1:7],
group = "2016",
fillOpacity =0.8,
color = ~pal2(argument[1:7]),
highlightOptions = highlightOptions(color = "black", weight = 2,
bringToFront = TRUE)) %>%
hideGroup("2016") %>%
#cliniccoordinates
addMarkers(
lng = cliniccoordinates$longitude, lat = cliniccoordinates$latitude,
label = cliniccoordinates$X,
labelOptions = labelOptions(noHide = F), group = cliniccoordinates$X) %>%
hideGroup(cliniccoordinates$X) %>%
addLegend("bottomleft", pal = pal2, values = argument,
title = argument2,
opacity = 1
) %>%
#Layers control
addLayersControl(
baseGroups = c("2004", "2006", "2009", "2012", "2016", "Nothing"),
overlayGroups = c(cliniccoordinates$X),
options = layersControlOptions(collapsed = TRUE)
)
}
colnames(arransimd)
[1] "DataZone" "LAName" "Rank" "Quintile" "Decile" "Vigintile" "Percentile" "IncRate"
[9] "IncRank" "EmpRate" "EmpRank" "HlthRank" "EduRank" "GAccRank" "HouseRank" "Shape_Leng"
[17] "Shape_Area" "year" "geometry" "lon" "lat" "listID"
functionmap(exampleshapes2$Rank, "Rank")
functionmap(exampleshapes2$Quintile, "Quintile")
functionmap(exampleshapes2$Quintile, "Decile")
functionmap(exampleshapes2$Vigintile, "Vigintile")
functionmap(exampleshapes2$Percentile, "Percentile")
functionmap(exampleshapes2$IncRate, "IncRate")
functionmap(exampleshapes2$IncRank, "IncRank")
functionmap(exampleshapes2$EmpRate, "EmpRate")
functionmap(exampleshapes2$HlthRank, "HlthRank")
functionmap(exampleshapes2$EduRank, "EduRank")
functionmap(exampleshapes2$GAccRank, "GAccRank")
functionmap(exampleshapes2$HouseRank, "HouseRank")